Prenatal phthalate exposure and cord blood DNA methylation

Exposure to phthalates has been shown to impede the human endocrine system, resulting in deleterious effects on pregnant women and their children. Phthalates modify DNA methylation patterns in infant cord blood. We examined the association between prenatal phthalate exposure and DNA methylation patterns in cord blood in a Korean birth cohort. Phthalate levels were measured in 274 maternal urine samples obtained during late pregnancy and 102 neonatal urine samples obtained at birth, and DNA methylation levels were measured in cord blood samples. For each infant in the cohort, associations between CpG methylation and both maternal and neonate phthalate levels were analyzed using linear mixed models. The results were combined with those from a meta-analysis of the levels of phthalates in maternal and neonatal urine samples, which were also analyzed for MEOHP, MEHHP, MnBP, and DEHP. This meta-analysis revealed significant associations between the methylation levels of CpG sites near the CHN2 and CUL3 genes, which were also associated with MEOHP and MnBP in neonatal urine. When the data were stratified by the sex of the infant, MnBP concentration was found to be associated with one CpG site near the OR2A2 and MEGF11 genes in female infants. In contrast, the concentrations of the three maternal phthalates showed no significant association with CpG site methylation. Furthermore, the data identified distinct differentially methylated regions in maternal and neonatal urine samples following exposure to phthalates. The CpGs with methylation levels that were positively associated with phthalate levels (particularly MEOHP and MnBP) were found to be enriched genes and related pathways. These results indicate that prenatal phthalate exposure is significantly associated with DNA methylation at multiple CpG sites. These alterations in DNA methylation may serve as biomarkers of maternal exposure to phthalates in infants and are potential candidates for investigating the mechanisms by which phthalates impact maternal and neonatal health.


Results
Subject characteristics. The subject characteristics of the mothers and the newborns for each center are presented in the following order: Ulsan, Ewha, and Dankook. The sample sizes for the maternal group were 159, 87, and 28, and those for the birth group were 74, 14, and 14, respectively. The average BMIs before pregnancy in the maternal group were 21.97 kg/m 2 (standard error (SE) = 0.25 kg/m 2 ), 22.42 kg/m 2 (SE = 0.34 kg/m 2 ), and 21.48 kg/m 2 (SE = 0.41 kg/m 2 ), respectively. In the birth group, the average BMIs were 22.38 kg/m 2 (SE = 0.36 kg/ m 2 ), 22.97 kg/m 2 (SE = 0.99 kg/m 2 ), and 21.39 kg/m 2 (SE = 0.73 kg/m 2 ), respectively. The number of current smokers in the maternal group was 3, 1, and 1, respectively, and there were no current smokers among the matched mothers for birth subjects in the birth group. In the birth group, the number of males was 43, 4, and 6, and for females, the numbers were 31, 10, and 8, respectively (Table 1) 35.13 µg/L, and 25.12 µg/L and in the birth group were 6.69 µg/L, 11.93 µg/L, and 9.16 µg/L. To test for differences in exposure levels among the three centers, Welch's one-way test 24 was performed on the log 2 -transformed data. For maternal group subjects, the test results for the concentrations of MEOHP, MEHHP, MnBP, and DEHP had p-values of 5.4E-03, 8.5E-03, 3.5E-01, and 5.1E-03, respectively, and for birth group subjects, the p-values were 8.3E-02, 1.3E-01, 2.5E-02, and 8.2E-02, respectively. Therefore, a meta-analysis was performed to combine the results from the three centers. The center-specific and total violin plots as well as the p-values of Welch's one-way test are shown in Fig. 1.

Detection of differentially methylated positions (DMPs).
After quality control, 793,925 CpGs were tested for their association with phthalate concentrations. Based on the FDR-adjusted 0.05 significance level using the Benjamini-Hochberg method and absolute value of weighted beta over 0.1, no CpGs were significantly associated with the four phthalate concentrations in the maternal group samples. Under the same threshold, cg12469381 and cg11229715 had a significant association with MEOHP and MnBP concentrations in the birth group samples, respectively. The center-specific and meta-combined regression coefficients, p-values, and annotation results are shown in Table 2. EWAS for DEHP exposure in the birth group also revealed one significant CpG, cg16419298, under the FDR-adjusted significance level of 0.05; however, the meta-combined regression coefficient was 0.024, which was not above the 0.1 threshold.
The sex-specific EWAS analysis revealed two significant CpGs associated with MnBP concentration in females in the birth group, cg05353481 and cg22290225, with an FDR-adjusted significance level of 0.05 (Table 2).
Differentially methylated regions (DMRs) for all samples. Regarding MEOHP, MEHHP, MnBP, and DEHP, 5, 5, 6, and 6 DMRs were identified, respectively, with 3, 5, 4, and 5 unique overlapping genes, respectively, for phthalate concentrations in the maternal group; similarly, for phthalate concentrations in the birth group, 13, 11, 5, and 13 DMRs were identified, respectively, with 11, 10, 2, and 11 unique overlapping genes, Gene set analysis (GSA). GSA was performed for MEOHP and MnBP concentrations in the birth group using their single significant DMPs. Under α = 0.05, for MEOHP and MnBP concentrations, two and 134 GO terms were identified, respectively, and two KEGG pathways were identified for MnBP concentration. The top ten terms are listed in Table 4. Under the FDR-adjusted significance level, no significant terms were identified. Previously identified DMRs were also used to perform gene set analyses. No terms were identified under the 0.05 FDR-adjusted significance level. Under α = 0.05, for MEOHP, MEHHP, MnBP, and DEHP concentrations in the maternal group, 22, 29, 5, and 29 GO terms were significant, respectively, and one KEGG pathway was identified for MEOHP concentration. Also under α = 0.05, for MEOHP, MEHHP, MnBP, and DEHP concentrationassociated terms in the birth group, 59, 41, 40, and 63 GO terms and 3, 0, 25, and 4 KEGG pathways, respectively, were identified (Supplementary Tables 2-5). The GO and KEGG terms associated with sex-specific DMRs are shown in Supplementary Tables 6-9 and 10-13 for females and males, respectively.

Discussion
In this study, we found evidence of the relationship between DNA methylation and prenatal phthalate exposure, during late pregnancy. Additionally, two CpG sites in female infants were significantly associated with urine phthalate concentrations at birth. We also identified a significant association between DMRs and maternal and neonatal urinary phthalate concentrations during late pregnancy with each phthalate type (MEHHP, MEOHP, DEHP, and MnBP). We found GO and KEGG pathways under α = 0.05, using the identified DMPs and DMRs associated with multiple genes enriched for pathways related to embryonic development and tumorigenesis, cell cycle progression, and signal transduction. Under the 0.05 FDR-adjusted significance level, no significant terms were identified. Individual p-values for each term are not independent, and if adjusted, the p-values can become smaller. Results can also be affected by the small sample sizes. Overall, our results suggest a sex-specific www.nature.com/scientificreports/ www.nature.com/scientificreports/ association between prenatal exposure to phthalates and DNA methylation at specific CpG sites, indicating that Korean infants may be differentially susceptible to phthalate-induced epigenetic alterations.
In the urine samples obtained after birth, MEOHP concentrations were associated with CpG sites (cg12469381) located near the chimerin 2 (CHN2) gene, which encodes GTPase-activating protein, found mainly in the pancreas and brain 25 . Another differentially methylated CpG site was identified at position 225,441,832 on chromosome 2, corresponding to cullin3 (CUL3), suggesting that CUL3 (cg11229715) also has a significant association with cord blood MnBP concentration. To the best of our knowledge, none of the genes associated with pregnancy phthalate concentrations covered in this study have been previously described in the context of epigenetic modifications or functions. However, based on the MRC-IEU EWAS catalog 26 , the CpG site cg12469381 extends beyond birth into late adolescence, impacting immune-neurodevelopmental functions, as shown in two population-based prospective birth cohorts 27 . Moreover, CUL3 is associated with pseudohypoaldosteronism and neurodevelopmental disorders, with or without autism or seizures 28 . Our results indicate that is a phthalate exposure-specific relationship with DNA methylation at specific CpG sites during pregnancy, but further research is needed to determine how prenatal exposure to phthalates impacts neurodevelopment and function.
In this study, the olfactory receptor family 2 subfamily A member 2 (OR2A2) gene for CpG sites (cg05353481) was significantly associated with birth urinary MnBP concentration, based on EWAS analysis for females. According to the EWAS catalog 26 , this gene is associated with age, maternal BMI, and nitrogen dioxide exposure 27,29,30 . The GO annotations associated with this gene include G protein-coupled receptor activity and taste receptor activity. Additionally, we found that in female infants another gene, cg22290225, annotated with Multiple EGF Like Domains 11 (MEGF11), was significantly associated with urine MnBP concentration. MEGF11 was manifested in the embryonic retina but was detected in horizontal cells and starburst amacrine cells in the first postnatal week and persisted into adulthood 31 . To date, none of the other gene studies have reported differential DNA methylation associated with urinary MnBP at birth. Sex-specific effects of some phthalates on methylation and expression of genes in cord blood have previously been reported 16,20 . However, the mechanisms involved in these sex-specific differences are still unclear. One plausible explanation is that the sex chromosome may contribute to sex-specific differences in metabolism. Previous studies have shown that the active X chromosome has more methylated CpGs in the body of genes 32,33 . Moreover, exposure to phthalates may negatively affect fetal development parameters, such as reduction in gestational age and preterm birth, and these effects can be sex-specific 4,13,34 . Considering this, it is important to note even small magnitudes of effects in these studies when they are sex-specific; our results provide further evidence for the association between prenatal phthalate exposure and epigenetic changes.
We identified more specific DMRs associated with phthalate metabolites in both maternal and infants, respectively, and found no overlap between them. Among the phthalate overlapping genes for MEOHP, MEHHP, and DEHP, in marternal, PLXNC1 belong to the family of trans-membrane receptors. It participates in the development of the nervous and immune systems through neuronal polarity, axon guidance, cellular motility, migration, and the immune response 35,36 . CHFR was significantly associated with maternal MnBP concentration in late pregnancy. As a tumorigenic factor, this gene promotes cell cycle progression and may induce ubiquitination by promoting polyubiquitination and degradation of HDAC1 37,38 . In addition, we examined the associations between DMR and each metabolite at birth. The EH domain-binding protein 1-like 1 (EHBP1L1) gene has a significantly hypomethylated DMR associated with MEOHP, MEHHP and DEHP concentrations at birth in infants. This gene is involved in the differentiation process of spermatogenesis and testicular activity 39 . A hypomethylated region of the microcephalin 1 (MCPH1) gene is involved in mitotic delay and irregular spindle orientation, resulting in defects in MCPH1 and an early change in cell division from symmetric to asymmetric, owing to the depletion of the progenitor pool and microencephaly 40,41 . Interestingly, urinary MnBP concentration was positively associated with two DMRs at birth. UVSSA encodes protein-coding genes and plays a role in the genome integrity homeostasis network. It is also involved in RNA polymerase II processing and repair factor recruitment. Defects in UVSSA affect the TC-NER pathway in RNAPII processing in various types of DNA damage, possibly contributing  www.nature.com/scientificreports/ to the development of a neurodegenerative phenotype common in disorders associated with genome instability 42 .
Overall, the DMRs found in our study in genes associated with phthalate exposure either in the urine of mothers during late pregnancy or infants are closely related to neurological functions, indicating the importance of the effects of phthalate exposure on genes associated with embryonic development. In our study, using gene set enrichment analysis, we showed that urinary MEOHP and MnBP concentrations at birth were associated with GO and KEGG pathways during late pregnancy. Several of these functional categories, such as regulation of signal transduction mediated by small GTPases and positive regulation of GTPase activity, were involved in biological processes associated with urinary MEOHP concentration at birth. In another study, a similar function was demonstrated for different phthalate metabolites, such as butyl benzyl phthalate (BBzP), which is associated with signal transduction mediated by small GTPases only in female infants 20 . In the present study, embryonic cleavage, positive regulation of mitotic metaphase, and mitotic sister chromatid   www.nature.com/scientificreports/ separation were significantly associated with urinary MnBP concentration at birth. A previous study suggested that phthalates may be associated with sperm methylation in or near genes important for early embryogenesis 43 . This would explain a biological pathway linking the observed inverse relationship between the antiandrogenic effects of phthalates and blastocyst quality. Additionally, we observed phthalate-induced methylation in sites associated with the Hedgehog signaling pathway or ubiquitin-mediated proteolysis. This signaling pathway has been demonstrated in vitro and in vivo system [43][44][45][46] . One study has reported that the effect of preconception exposure to DEHP on genome-wide DNA methylation and gene expression profiles in mice, was related to the Hedgehog signaling pathway 43 . In addition, other studies have found that in utero exposure to dibutyl phthalate (DBP) led to abnormal proliferation of testicular Sertoli cells in prepubertal mice by modulating the ubiquitination of the key proliferation-related protein IRAK1 via the downregulation of Peli2 46 . Several studies have shown that phthalate exposure was associated with numerous enriched inflammatory pathways (e.g., NF-κB, MAPK, TNFβ) 19,47 . In our study, these pathways were not significantly associated with phthalate exposure in late pregnancy or phthalate levels in cord blood. However, methylation is positively associated with the expression of some genes, and it is necessary to investigate whether the methylation process is involved in the regulation of these pathways. Understanding the molecular differences affecting these CpGs sites may lead to a better understanding of possible abnormalities in fetal development caused by prenatal phthalate exposure.
Our study was specifically designed to investigate the relationship between prenatal phthalate exposure and DNA methylation in cord blood. However, this study has some limitations. First, the sample size of cord blood was relatively small, which limits the statistical power of the analysis or the availability of methods. The metaanalysis revealed two CpG sites and two other CpG sites in females that were significantly associated with cord blood concentration. However, we could not find any correlation between maternal phthalate concentrations and CpG sites. A larger sample size would allow for greater statistical power and could help identify relevant DMPs. Further studies with a larger population are needed to confirm the findings of our study. Second, the participants in this study were in the late stages of pregnancy. Phthalates have an extremely short half-life and are rapidly excreted from the body 48 ; therefore, a single estimate is not representative of the entire pregnancy. However, most significant DMRs are generally observed in late pregnancy 14,20 . Moreover, we defined the possible biological processes affected by prenatal phthalate exposure based on GO and KEGG pathways, which may be useful for further research. Lastly, a previous study suggested that some probes in the Illumina MethylationEPIC BeadChip can bind with polymorphic sites or hybridize with non-specific sites. While these types of probes can still represent methylation status at the appropriate sites, the results could be inaccurate if the probe itself is methylated 49 . To better explore the epigenome associated with phthalates, our study used the EPIC BeadChip (loci:850,000) instead of the 450 k BeadChips (loci:450,000). Although the EPIC BeadChip approach has become the standard in EWAS studies, the reported percentages of differentially methylated CpGs and enrichment-associated genes may be biased because of the lower coverage of epigenome sequencing. Therefore, these probes need to be interpreted more carefully, and additional reporting or biological validation in different cohorts may be required.
In conclusion, there is a significant association between prenatal phthalate exposure and DNA methylation at multiple CpG sites. Additionally, urinary MnBP levels at birth were significantly associated with CpG sites in female infants. Our study has uncovered several functional mechanisms and associated genes, providing insight into the effects of phthalate exposure on fetal development and neurodevelopmental disorders. Further research is needed to investigate the mechanisms behind the adverse effects of prenatal phthalate exposure on human well-being.

Methods
Study population and phthalate measurement. From 2006 to 2010, pregnant women within 20 weeks of pregnancy and over 18 years of age were recruited from the Mother and Child's Environmental Health (MOCEH) multicenter prospective cohort study. The health centers were located in Ulsan, Ewha, and Dankook. Patients who were planning to move out of the locations within a year or were diagnosed with cognitive impairment or mental disease were excluded from the study. Details on the recruitment process for this cohort can be found in a previous study 50 . All participants provided written informed consent during their first prenatal visits, and during routine face-to-face visits. Self-administered questionnaires were completed to collect sociodemographic information. This study was approved by the Institutional Review Board of Kangwon National University Hospital (2017-11-006). This study was conformed to the tenets of the Declaration of Helsinki.
Phthalate concentrations in the urine samples were measured according to the MOCEH study protocol 50 . Samples were analyzed using high-performance liquid chromatography-tandem mass spectrometry (Agilent Technologies, Inc., California, USA). The limits of detection (LODs) for MEOHP, MEHHP, and MnBP were 0.56 µg/L, 0.49 µg/L, and 0.44 µg/L, respectively. Measurements below the detection limit were recorded as half the LOD. The DEHP concentration was calculated by summing MEOHP and MEHHP concentrations 51 . DNA methylation data preprocessing. DNA  www.nature.com/scientificreports/ essed following the recommended quality control (QC) workflow of the 'ewastools' package in R 54,55 . Probes with detection p-values above the significance threshold of 0.01 were set to missing, and dye-bias was corrected with the RELIC method 56 using the Theil-Sen estimator. Then, beta values were calculated as the ratio of methylated to total (methylated and unmethylated combined) signal intensities. Leukocyte composition was estimated using the beta values with the Houseman method 57 . Based on the reference panel given by Salas et al., cell deconvolution was achieved using the CpGs optimized with Identifying Optimal Libraries (IDOL) and the "estimateCellCounts2" function with the "CordBlood" option in the "FlowSorted.Blood.EPIC" package in R, which is an expansion of the "minfi" package 58-60 . Quality control (QC). The "ewastools" package was used to perform QC on the subjects and CpGs. For the subject-level QC, following the methods of Park et al 53 ., subjects were excluded according to the following criteria: (1) one of the twins; (2) recommended thresholds of the 17 metrics from the Illumina BeadArray Controls Reporter Software Guide 61 ; (3) inconsistency between observed sex and sex inferred from normalized X/Y probe intensities; (4) outliers or duplicates for the SNP probe intensities; (5) outliers for the principal components of the autosomal beta values; or (6) outliers for the leukocyte composition. For step (4), a detailed quality check procedure 54 was performed according to a previously described process. The beta values for SNP probes were trained with a mixture model with three beta distributions for subjects with genotypes AA/Aa/aa and a uniform distribution for outlying subjects. The average log odd of being outliers was calculated for each subject across all SNPs, and if this value was above -4, the corresponding subject was marked as an outlier. Moreover, the posterior probabilities for each genotype were used to calculate the proportions of SNPs for which both had the same genotype for each pair of subjects, excluding those with proportions exceeding 90%. For steps (5) and (6), a multidimensional scaling (MDS) plot was generated, and outliers were manually detected. Likewise, leukocyte compositions were plotted separately for each cell type, and manually-detected outliers were filtered out. The number of subjects excluded by QC procedures is shown in Supplementary Covariate imputation. After excluding one randomly selected twin, 383 subjects were used for the imputation of missing values of the following seven covariates using the "missForest" package in R 55,63 : infant sex, maternal pre-pregnancy body mass index (BMI), smoking history, depending on whether the participants had smoked at least 400 cigarettes, and early and late pregnancy cotinine and creatinine concentrations in maternal blood. For the dataset used for the analysis, please refer to Supplementary Table 2 of a previous published study 55 .
Using imputed values for smoking history and late pregnancy cotinine concentrations, participants were classified as "current smokers" if their late pregnancy cotinine concentrations (measured or imputed) were not less than 200 µg/g or if they answered "yes" to a questionnaire on whether they continued to smoke. For the actual statistical analysis, subjects with missing sex information were excluded because (1) the QC steps require sex to be reported and (2) it would be more robust to use only observed values if stratification by sex was performed.

Epigenome-wide association study (EWAS).
An EWAS was performed on MEOHP, MEHHP, MnBP, and DEHP concentrations in the maternal and birth groups. EWAS was initially performed separately for the three centers, and the results were later combined with the meta-analysis in order to minimize the heterogeneity among centers. For each center, linear regression was performed on each of the CpG sites using the "limma" package in R to test the association between the beta value of methylation level and log2-transformed phthalate concentration in urine 55,64 . The covariates included in the model were the study enrollment year, infant sex, maternal pre-pregnancy BMI, current smoking status, average monthly household income, and estimated leukocyte composition. Batch information was included as a random effect to adjust for batch effects. The p-values from each center were combined using Liptak's method 65 , and regression coefficients were combined using weights inversely proportional to the standard errors. Also, using the sex of birth subjects, sex-stratified EWAS was performed for samples in the birth group. Due to the limited sample size, samples from the three centers were combined and analyzed together (mega-analysis). The sex-specific sample sizes for each center are listed in Table 1.
Differentially methylated region (DMR) analysis. The DMRs associated with phthalate concentration in urine were separately identified from the three centers and combined with meta-analysis using the "dmrff " package in R and Rex (version 3.6.0) 66,67 . The methylation levels of each CpG site were first transformed using inverse normal transformation before the analysis to obtain results robust to outliers and the normal distribution assumption. Regions with a maximum distance of 500 bp between consecutive features and at least two significant probes at a significance level of 0.05 were identified. The identified DMRs were evaluated using the 0.05 Bonferroni adjusted significance level. ENSEMBL_MART_ENSEMBL BioMart database and the hsapiens_ gene_ensembl database in the Ensembl genome browser (version: GRCh37) were used for annotation. Megaanalysis for sex-specific DMR identification in the birth group was performed using the "dmrff " R package. Annotation was performed using the same procedure as described above.

Gene set analysis (GSA). Gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes
(KEGG) pathways were identified using the significant CpGs and DMRs. Terms with at least five CpG sites were used to create a gene set, and each gene set was tested using the 0.05 FDR-adjusted significance level. GSA was